**August 26, 2019
**Code corrections to PDKM's do files by D.M. Gibler
**Note: The code that follows is the code used by PDKM to generate replication of the Braithwaite and Lemke piece.
**Any changes we made to the code are prefaced by "***dmg:..***".  All file calls were of course changed
**because we use our own directories. Comments that begin with only one asterisk (*) are in the original code.


* This file is a re-replication of Gibler, Miller, and Little replication of Braithwaite and Lemke (2011). 
* This is the second of two Braithwaite and Lemke (BL) replication files.
* It focuses on replicating Braithwaite and Lemke's statistical analysis.

clear
capture log close
log using "Braithwaite and Lemke Analysis Replication.log", replace
set more off

* Merge in variables created in the previous do file.

*** dmg: (cut) use "B&L_CMPS.dta", clear
* Note: This is the BL replication dataset from Gibler's dataverse, but we added the variables ccodea and ccodeb and dropped some estimated variables.

*** dmg: The <base data> is the corrected B&L data described in the other do file.  
use "base data.dta", clear

*** dmg: We create ccodea and ccodeb in the reconstructed data directly and comment out PDKM code.
gen ccodea=ccode1
gen ccodeb=ccode2

* The dyadid is missing in some cases, so we fix this first.
gen ccodehigh=ccodea if ccodea>ccodeb
replace ccodehigh=ccodeb if ccodeb>ccodea
gen ccodelow=ccodea if ccodea<ccodeb
replace ccodelow=ccodeb if ccodeb<ccodea
***dmg: (cut) replace dyadid=ccodelow*1000+ccodehigh
gen dyadid=ccodelow*1000+ccodehigh

merge m:1 dyadid year using "BL Dyadic MIDs 3.0.dta"
drop if _merge==2
drop _merge
replace onset_mid3=0 if onset_mid3==.

merge m:1 dyadid year using "BL Dyadic MIDs 4.3.dta"
drop if _merge==2
drop _merge
replace onset_mid4=0 if onset_mid4==.

merge m:1 dyadid year using "BL Dyadic MIDs GML ISQ.dta"
drop if _merge==2
drop _merge
replace onset_gml1=0 if onset_gml1==.

*** dmg: (cut) merge m:1 dyadid year using "BL Dyadic MIDs GML 2.1.dta"
*** dmg: (cut) drop if _merge==2
*** dmg: (cut) drop _merge
*** dmg: (cut) replace onset_gml2=0 if onset_gml2==.

*** dmg: The following renames variables in the reconstructed data to fit PDKM's code.
rename riv1 thom_rival
rename bothminors min_min 
rename democracy joint_democ
rename satisfaction jointsatis
rename hadispute terrdisp
rename parity powerratio
rename contiguity contiguous	

*** dmg: The base data included only the independent variables, no dispute variables.  Therefore, we have
*** dmg: to label the territorial disputes.  The following do file is a simple script that lavels all
*** dmg: territorial disputes.
do "territorial dispute labeler.do"

* Labeling
label var thom_rival "Rivalry"
label var min_min "Both Minors"
label var joint_democ "Joint Democracy"
label var jointsatis "Joint Satisfaction"
label var terrdisp "Territorial Claim"
label var powerratio "Power Preponderance"
label var defense "Allied" 
label var territory "Territorial Dispute"
label var contiguous "Contiguous"				

*** dmg: We reconstructed the data as directed-dyad data.  B&L uses non-directed data, so we drop directionality.
drop if ccode1>ccode2

*** dmg: PDKM just used the original peace years for temporal duration adjustment in every single model even though
*** dmg: the dependent variables changed substantially across datasets.
*** dmg: We correct this by generatin separate peaceyears and splines using the ado BTSCS
btscs onset_mid3 year ccode1 ccode2, g(peaceyearsmid3)
gen peaceyearsmid3sqr=peaceyearsmid3^2
gen peaceyearsmid3cub=peaceyearsmid3^3

btscs onset_gml1 year ccode1 ccode2, g(peaceyearsgml)
gen peaceyearsgmlsqr=peaceyearsgml^2
gen peaceyearsgmlcub=peaceyearsgml^3

btscs onset_mid4 year ccode1 ccode2, g(peaceyears43)
gen peaceyears43sqr=peaceyears43^2
gen peaceyears43cub=peaceyears43^3

****************************************  ANALYSIS  **************************************************************************

* ANALYSIS WITH MID 3.0 DATA (ORIGINALLY USED BY BL)

* In addition to running the same regressions as BL, we add some matrix notation to help us make the figures.

* MID 3.0 Recip Model

*** dmg: (cut) heckprob recip_mid3 joint_democ territory jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob recip_mid3 joint_democ territory jointsatis powerratio defense , ///
	sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyearsmid3*) cluster(dyadid)
	
eststo m1_bl
***
matrix m1_bl_out = J(5,2,.)
matrix coln m1_bl_out = coef se
matrix rown m1_bl_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m1_bl_out[`i',1] = beta[1,`i']
	matrix m1_bl_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m1_bl_sel = J(9,2,.)
matrix coln m1_bl_sel = coef se
matrix rown m1_bl_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m1_bl_sel[`i'-6,1] = beta[1,`i']
	matrix m1_bl_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m1_bl_sel[9,1] = beta[1,19]
matrix m1_bl_sel[9,2] = sqrt(cov[19,19])
***

* MID 3.0 Use of Force Model

*** dmg: (cut) heckprob force_mid3 joint_democ territory  jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob force_mid3 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyearsmid3*) cluster(dyadid)
	
eststo m2_bl
***
matrix m2_bl_out = J(5,2,.)
matrix coln m2_bl_out = coef se
matrix rown m2_bl_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m2_bl_out[`i',1] = beta[1,`i']
	matrix m2_bl_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m2_bl_sel = J(9,2,.)
matrix coln m2_bl_sel = coef se
matrix rown m2_bl_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m2_bl_sel[`i'-6,1] = beta[1,`i']
	matrix m2_bl_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m2_bl_sel[9,1] = beta[1,19]
matrix m2_bl_sel[9,2] = sqrt(cov[19,19])
***


* MID 3.0 Mutual Use of Force Model

*** dmg: (cut) heckprob mut_force_mid3 joint_democ territory  jointsatis powerratio defense  ///
*** dmg: (cut) 		,sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob mut_force_mid3 joint_democ territory  jointsatis powerratio defense  ///
	,sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyearsmid3*) cluster(dyadid)
	
eststo m3_bl
***
matrix m3_bl_out = J(5,2,.)
matrix coln m3_bl_out = coef se
matrix rown m3_bl_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m3_bl_out[`i',1] = beta[1,`i']
	matrix m3_bl_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m3_bl_sel = J(9,2,.)
matrix coln m3_bl_sel = coef se
matrix rown m3_bl_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m3_bl_sel[`i'-6,1] = beta[1,`i']
	matrix m3_bl_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m3_bl_sel[9,1] = beta[1,19]
matrix m3_bl_sel[9,2] = sqrt(cov[19,19])
***


* MID 3.0 Fatal Model

*** dmg: (cut) heckprob fatal_mid3 joint_democ territory  jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob fatal_mid3 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyearsmid3*) cluster(dyadid)
	
eststo m4_bl
***
matrix m4_bl_out = J(5,2,.)
matrix coln m4_bl_out = coef se
matrix rown m4_bl_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m4_bl_out[`i',1] = beta[1,`i']
	matrix m4_bl_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m4_bl_sel = J(9,2,.)
matrix coln m4_bl_sel = coef se
matrix rown m4_bl_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m4_bl_sel[`i'-6,1] = beta[1,`i']
	matrix m4_bl_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m4_bl_sel[9,1] = beta[1,19]
matrix m4_bl_sel[9,2] = sqrt(cov[19,19])
***


* MID 3.0 High Fatality Model

*** dmg: (cut) heckprob high_fatal_mid3 joint_democ territory  jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis  ///
*** dmg: (cut) 		terrdisp  powerratio defense peaceyears pax2 pax3) cluster(dyadid)

heckprob high_fatal_mid3 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis  ///
	terrdisp  powerratio defense peaceyearsmid3*) cluster(dyadid)
	
eststo m5_bl
***
matrix m5_bl_out = J(5,2,.)
matrix coln m5_bl_out = coef se
matrix rown m5_bl_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m5_bl_out[`i',1] = beta[1,`i']
	matrix m5_bl_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m5_bl_sel = J(9,2,.)
matrix coln m5_bl_sel = coef se
matrix rown m5_bl_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m5_bl_sel[`i'-6,1] = beta[1,`i']
	matrix m5_bl_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m5_bl_sel[9,1] = beta[1,19]
matrix m5_bl_sel[9,2] = sqrt(cov[19,19])
***


* MID 3.0 War Model

*** dmg: (cut) heckprob war_mid3 joint_democ territory  jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)

heckprob war_mid3 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid3 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyearsmid3*) cluster(dyadid)
	
eststo m6_bl
***
matrix m6_bl_out = J(5,2,.)
matrix coln m6_bl_out = coef se
matrix rown m6_bl_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m6_bl_out[`i',1] = beta[1,`i']
	matrix m6_bl_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m6_bl_sel = J(9,2,.)
matrix coln m6_bl_sel = coef se
matrix rown m6_bl_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m6_bl_sel[`i'-6,1] = beta[1,`i']
	matrix m6_bl_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m6_bl_sel[9,1] = beta[1,19]
matrix m6_bl_sel[9,2] = sqrt(cov[19,19])
***


*************************************************************************************************************************************************

* ANALYSIS WITH GML 2.1 DATA


* We use version 2.1 of the GML data to be consistent with our analysis of Weeks (for which using version 2.1 was necessary to obtain incident-level data).
* Results using the GML V1.7 dataset that was used in the ISQ piece can be seen by substituting “gml1” for “gml2” in the regressions below.

*** dmg: We substitute gml1 below for ease of code, but they're now equivalent given last do file.

* GML 2.1 Recip Model

*** dmg: (cut) heckprob recip_gml2 joint_democ territory jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_gml2 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob recip_gml1 joint_democ territory jointsatis powerratio defense , ///
	sel(onset_gml1 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyearsgml*) cluster(dyadid)
eststo m1_bl
***
matrix m1_gml_out = J(5,2,.)
matrix coln m1_gml_out = coef se
matrix rown m1_gml_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m1_gml_out[`i',1] = beta[1,`i']
	matrix m1_gml_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m1_gml_sel = J(9,2,.)
matrix coln m1_gml_sel = coef se
matrix rown m1_gml_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m1_gml_sel[`i'-6,1] = beta[1,`i']
	matrix m1_gml_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m1_gml_sel[9,1] = beta[1,19]
matrix m1_gml_sel[9,2] = sqrt(cov[19,19])
***

* GML 2.1 Use of Force Model

*** dmg: (cut) heckprob force_gml2 joint_democ territory  jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_gml2 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob force_gml1 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_gml1 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyearsgml*) cluster(dyadid)
	
eststo m2_bl
***
matrix m2_gml_out = J(5,2,.)
matrix coln m2_gml_out = coef se
matrix rown m2_gml_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m2_gml_out[`i',1] = beta[1,`i']
	matrix m2_gml_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m2_gml_sel = J(9,2,.)
matrix coln m2_gml_sel = coef se
matrix rown m2_gml_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m2_gml_sel[`i'-6,1] = beta[1,`i']
	matrix m2_gml_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m2_gml_sel[9,1] = beta[1,19]
matrix m2_gml_sel[9,2] = sqrt(cov[19,19])
***


* GML 2.1 Mutual Use of Force Model

*** dmg: (cut) heckprob mut_force_gml2 joint_democ territory  jointsatis powerratio defense  ///
*** dmg: (cut) 		,sel(onset_gml2 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)

heckprob mut_force_gml1 joint_democ territory  jointsatis powerratio defense  ///
	,sel(onset_gml1 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyearsgml*) cluster(dyadid)
	
eststo m3_bl
***
matrix m3_gml_out = J(5,2,.)
matrix coln m3_gml_out = coef se
matrix rown m3_gml_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m3_gml_out[`i',1] = beta[1,`i']
	matrix m3_gml_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m3_gml_sel = J(9,2,.)
matrix coln m3_gml_sel = coef se
matrix rown m3_gml_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m3_gml_sel[`i'-6,1] = beta[1,`i']
	matrix m3_gml_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m3_gml_sel[9,1] = beta[1,19]
matrix m3_gml_sel[9,2] = sqrt(cov[19,19])
***


* GML 2.1 Fatal Model

*** dmg: (cut) heckprob fatal_gml2 joint_democ territory  jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_gml2 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob fatal_gml1 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_gml1 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyearsgml*) cluster(dyadid)
eststo m4_bl
***
matrix m4_gml_out = J(5,2,.)
matrix coln m4_gml_out = coef se
matrix rown m4_gml_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m4_gml_out[`i',1] = beta[1,`i']
	matrix m4_gml_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m4_gml_sel = J(9,2,.)
matrix coln m4_gml_sel = coef se
matrix rown m4_gml_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m4_gml_sel[`i'-6,1] = beta[1,`i']
	matrix m4_gml_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m4_gml_sel[9,1] = beta[1,19]
matrix m4_gml_sel[9,2] = sqrt(cov[19,19])
***


* GML 2.1 High Fatality Model
*** dmg: (cut) heckprob high_fatal_gml2 joint_democ territory  jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_gml2 = contiguous thom_rival min_min joint_democ  jointsatis  ///
*** dmg: (cut) 		terrdisp  powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob high_fatal_gml1 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_gml1 = contiguous thom_rival min_min joint_democ  jointsatis  ///
	terrdisp  powerratio defense peaceyearsgml*) cluster(dyadid)
eststo m5_bl
***
matrix m5_gml_out = J(5,2,.)
matrix coln m5_gml_out = coef se
matrix rown m5_gml_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m5_gml_out[`i',1] = beta[1,`i']
	matrix m5_gml_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m5_gml_sel = J(9,2,.)
matrix coln m5_gml_sel = coef se
matrix rown m5_gml_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m5_gml_sel[`i'-6,1] = beta[1,`i']
	matrix m5_gml_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m5_gml_sel[9,1] = beta[1,19]
matrix m5_gml_sel[9,2] = sqrt(cov[19,19])
***


* GML 2.1 War Model
*** dmg: (cut) heckprob war_gml2 joint_democ territory  jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_gml2 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob war_gml1 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_gml1 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyearsgml*) cluster(dyadid)
eststo m6_bl
***
matrix m6_gml_out = J(5,2,.)
matrix coln m6_gml_out = coef se
matrix rown m6_gml_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m6_gml_out[`i',1] = beta[1,`i']
	matrix m6_gml_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m6_gml_sel = J(9,2,.)
matrix coln m6_gml_sel = coef se
matrix rown m6_gml_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m6_gml_sel[`i'-6,1] = beta[1,`i']
	matrix m6_gml_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m6_gml_sel[9,1] = beta[1,19]
matrix m6_gml_sel[9,2] = sqrt(cov[19,19])
***



*************************************************************************************************************************************************

* ANALYSIS WITH MID 4.3 DATA

* MID 4.3 Recip Model

*** dmg: (cut) heckprob recip_mid4 joint_democ territory jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob recip_mid4 joint_democ territory jointsatis powerratio defense , ///
	sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears43*) cluster(dyadid)
eststo m1_bl
***
matrix m1_mid4_out = J(5,2,.)
matrix coln m1_mid4_out = coef se
matrix rown m1_mid4_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m1_mid4_out[`i',1] = beta[1,`i']
	matrix m1_mid4_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m1_mid4_sel = J(9,2,.)
matrix coln m1_mid4_sel = coef se
matrix rown m1_mid4_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m1_mid4_sel[`i'-6,1] = beta[1,`i']
	matrix m1_mid4_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m1_mid4_sel[9,1] = beta[1,19]
matrix m1_mid4_sel[9,2] = sqrt(cov[19,19])
***

* MID 4.3 Use of Force Model

*** dmg: (cut) heckprob force_mid4 joint_democ territory  jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob force_mid4 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears43*) cluster(dyadid)
eststo m2_bl
***
matrix m2_mid4_out = J(5,2,.)
matrix coln m2_mid4_out = coef se
matrix rown m2_mid4_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m2_mid4_out[`i',1] = beta[1,`i']
	matrix m2_mid4_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m2_mid4_sel = J(9,2,.)
matrix coln m2_mid4_sel = coef se
matrix rown m2_mid4_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m2_mid4_sel[`i'-6,1] = beta[1,`i']
	matrix m2_mid4_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m2_mid4_sel[9,1] = beta[1,19]
matrix m2_mid4_sel[9,2] = sqrt(cov[19,19])
***


* MID 4.3 Mutual Use of Force Model

*** dmg: (cut) heckprob mut_force_mid4 joint_democ territory  jointsatis powerratio defense  ///
*** dmg: (cut) 		,sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob mut_force_mid4 joint_democ territory  jointsatis powerratio defense  ///
	,sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears43*) cluster(dyadid)
eststo m3_bl
***
matrix m3_mid4_out = J(5,2,.)
matrix coln m3_mid4_out = coef se
matrix rown m3_mid4_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m3_mid4_out[`i',1] = beta[1,`i']
	matrix m3_mid4_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m3_mid4_sel = J(9,2,.)
matrix coln m3_mid4_sel = coef se
matrix rown m3_mid4_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m3_mid4_sel[`i'-6,1] = beta[1,`i']
	matrix m3_mid4_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m3_mid4_sel[9,1] = beta[1,19]
matrix m3_mid4_sel[9,2] = sqrt(cov[19,19])
***


* MID 4.3 Fatal Model

*** dmg: (cut) heckprob fatal_mid4 joint_democ territory  jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob fatal_mid4 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears43*) cluster(dyadid)
eststo m4_bl
***
matrix m4_mid4_out = J(5,2,.)
matrix coln m4_mid4_out = coef se
matrix rown m4_mid4_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m4_mid4_out[`i',1] = beta[1,`i']
	matrix m4_mid4_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m4_mid4_sel = J(9,2,.)
matrix coln m4_mid4_sel = coef se
matrix rown m4_mid4_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m4_mid4_sel[`i'-6,1] = beta[1,`i']
	matrix m4_mid4_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m4_mid4_sel[9,1] = beta[1,19]
matrix m4_mid4_sel[9,2] = sqrt(cov[19,19])
***


* MID 4.3 High Fatality Model

*** dmg: (cut) heckprob high_fatal_mid4 joint_democ territory  jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis  ///
*** dmg: (cut) 		terrdisp  powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob high_fatal_mid4 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis  ///
	terrdisp  powerratio defense peaceyears43*) cluster(dyadid)
eststo m5_bl
***
matrix m5_mid4_out = J(5,2,.)
matrix coln m5_mid4_out = coef se
matrix rown m5_mid4_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m5_mid4_out[`i',1] = beta[1,`i']
	matrix m5_mid4_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m5_mid4_sel = J(9,2,.)
matrix coln m5_mid4_sel = coef se
matrix rown m5_mid4_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m5_mid4_sel[`i'-6,1] = beta[1,`i']
	matrix m5_mid4_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m5_mid4_sel[9,1] = beta[1,19]
matrix m5_mid4_sel[9,2] = sqrt(cov[19,19])
***


* MID 4.3 War Model

*** dmg: (cut) heckprob war_mid4 joint_democ territory  jointsatis powerratio defense , ///
*** dmg: (cut) 		sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
*** dmg: (cut) 		terrdisp powerratio defense peaceyears pax2 pax3) cluster(dyadid)
heckprob war_mid4 joint_democ territory  jointsatis powerratio defense , ///
	sel(onset_mid4 = contiguous thom_rival min_min joint_democ  jointsatis ///
	terrdisp powerratio defense peaceyears43*) cluster(dyadid)
eststo m6_bl
***
matrix m6_mid4_out = J(5,2,.)
matrix coln m6_mid4_out = coef se
matrix rown m6_mid4_out = joint_democ territory jointsatis powerratio defense
matrix beta = e(b)
matrix cov = e(V)

foreach i of num 1/5 {
	matrix m6_mid4_out[`i',1] = beta[1,`i']
	matrix m6_mid4_out[`i',2] = sqrt(cov[`i',`i'])
}

matrix m6_mid4_sel = J(9,2,.)
matrix coln m6_mid4_sel = coef se
matrix rown m6_mid4_sel = contiguous thom_rival min_min joint_democ jointsatis ///
						terrdisp powerratio defense Rho
*Pull each variable, plus Rho
foreach i of num 7/14 {
	matrix m6_mid4_sel[`i'-6,1] = beta[1,`i']
	matrix m6_mid4_sel[`i'-6,2] = sqrt(cov[`i',`i'])
}
matrix m6_mid4_sel[9,1] = beta[1,19]
matrix m6_mid4_sel[9,2] = sqrt(cov[19,19])
***

*******************************************************************************************************************************

* MAKE FIGURES

coefplot ///
	(matrix(m1_bl_sel[.,1]), se(m1_mid4_sel[.,2]) m(smcircle) mlw(medium) msize(medlarge) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m1_gml_sel[.,1]), se(m1_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m1_mid4_sel[.,1]), se(m1_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.2)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(col(3) ring()  region(lcolor(none)) size(2)) ///
	xtitle("Coefficient Value" , size(small)) ///
	ytitle(" ", size(small) margin(b=0 l=0 t=0 r=0) ) ///
	subtitle("Onset" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(-1 -.5 0 .5 1 , labsize(small) angle(0)  ) ///
	ylabel(,notick   labsize(small)  ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(20) xsize() fysize() ///
		graphregion(margin(l=0 r=-2)) ///
	levels(95 90) saving(m1_sel, replace) 
coefplot ///
	(matrix(m1_bl_out[.,1]), se(m1_bl_out[.,2]) m(smcircle) mlw(medium) msize(medium) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m1_gml_out[.,1]), se(m1_gml_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m1_mid4_out[.,1]), se(m1_mid4_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.2)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(off) ///
	xtitle("") ///
	ytitle(" ", size(small) margin(b=0 l=5.5 t=0 r=0)) ///
	title("Model 1"  , size(medsmall)) ///
	subtitle("Reciprocation" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noline noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(none  ) ///
	ylabel(,notick   labsize(small)  ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(5) xsize() fysize(80) ///
	graphregion(margin(l=0 r=-2 b=-8)) ///
	levels(95 90) saving(m1_out, replace)
graph combine m1_out.gph m1_sel.gph,  xsize(3) ysize(6) iscale(1.2) ///
	 scheme(s1mono) rows(2) saving(m1, replace) 

	 
	 
coefplot ///
	(matrix(m2_bl_sel[.,1]), se(m2_mid4_sel[.,2]) m(smcircle) mlw(medium) msize(medlarge) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m2_gml_sel[.,1]), se(m2_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m2_mid4_sel[.,1]), se(m2_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.2)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(col(3) ring()  region(lcolor(none)) size(2)) ///
	xtitle("Coefficient Value" , size(small)) ///
	ytitle(" ", size(small) margin(0 15 0 0) ) ///
	subtitle("Onset" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(-1 -.5 0 .5 1 , labsize(small) angle(0)  ) ///
	ylabel(none ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(20) xsize() fysize() ///
	graphregion(margin(l=0 r=18)) ///
	levels(95 90) saving(m2_sel, replace) 
coefplot ///
	(matrix(m2_bl_out[.,1]), se(m2_bl_out[.,2]) m(smcircle) mlw(medium) msize(medium) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m2_gml_out[.,1]), se(m2_gml_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m2_mid4_out[.,1]), se(m2_mid4_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.2)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(off) ///
	xtitle("") ///
	ytitle(" ", size(small) margin(0 15 0 0)) ///
	title("Model 2"  , size(medsmall)) ///
	subtitle("Use of Force" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noline noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(none) ///
	ylabel(none) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(5) xsize() fysize(80) ///
	graphregion(margin(l=0 r=18 b=-8)) ///
	levels(95 90) saving(m2_out, replace)

graph combine m2_out.gph m2_sel.gph,  xsize(3) ysize(6) iscale(1.2) ///
	 scheme(s1mono) rows(2) saving(m2, replace) 
	 
	 
*Model 3
coefplot ///
	(matrix(m3_bl_sel[.,1]), se(m3_mid4_sel[.,2]) m(smcircle) mlw(medium) msize(medlarge) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m3_gml_sel[.,1]), se(m3_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m3_mid4_sel[.,1]), se(m3_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.2)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(col(3) ring()  region(lcolor(none)) size(2)) ///
	xtitle("Coefficient Value" , size(small)) ///
	ytitle(" ", size(small) ) ///
	subtitle("Onset" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(-1 -.5 0 .5 1 , labsize(small) angle(0)  ) ///
	ylabel(none ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(20) xsize() fysize() ///
	graphregion(margin(l=9 r=23.5)) ///
	levels(95 90) saving(m3_sel, replace) 
coefplot ///
	(matrix(m3_bl_out[.,1]), se(m3_bl_out[.,2]) m(smcircle) mlw(medium) msize(medium) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m3_gml_out[.,1]), se(m3_gml_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m3_mid4_out[.,1]), se(m3_mid4_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.2)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(off) ///
	xtitle("") ///
	ytitle(" ", size(small) margin(0 0 0 0)) ///
	title("Model 3"  , size(medsmall)) ///
	subtitle("Mutual Use of Force" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noline noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(none  ) ///
	ylabel(none ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(5) xsize() fysize(80) ///
	graphregion(margin(l=9 r=23.5 b=-8)) ///
	levels(95 90) saving(m3_out, replace)
	 
graph combine m3_out.gph m3_sel.gph,  xsize(3) ysize(6) iscale(1.2) ///
	 scheme(s1mono) rows(2) saving(m3, replace) 
	

*Model4
coefplot ///
	(matrix(m4_bl_sel[.,1]), se(m4_mid4_sel[.,2]) m(smcircle) mlw(medium) msize(medlarge) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m4_gml_sel[.,1]), se(m4_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m4_mid4_sel[.,1]), se(m4_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.2)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(col(3) ring()  region(lcolor(none)) size(2)) ///
	xtitle("Coefficient Value" , size(small)) ///
	ytitle(" ", size(small) ) ///
	subtitle("Onset" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(-1 -.5 0 .5 1 , labsize(small) angle(0)  ) ///
	ylabel(,notick   labsize(small)  ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(20) xsize() fysize() ///
	graphregion(margin(l=0 r=-2)) ///
	levels(95 90) saving(m4_sel, replace) 
coefplot ///
	(matrix(m4_bl_out[.,1]), se(m4_bl_out[.,2]) m(smcircle) mlw(medium) msize(medium) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m4_gml_out[.,1]), se(m4_gml_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m4_mid4_out[.,1]), se(m4_mid4_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.2)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(off) ///
	xtitle("") ///
	ytitle(" ", size(small) margin(0 5.5 0 0)) ///
	title("Model 4"  , size(medsmall)) ///
	subtitle("Fatal MID" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noline noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(none  ) ///
	ylabel(,notick   labsize(small)  ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(5) xsize() fysize(80) ///
	graphregion(margin(l=0 r=-2 b=-8)) ///
	levels(95 90) saving(m4_out, replace)
graph combine m4_out.gph m4_sel.gph,  xsize(3) ysize(6) iscale(1.2) ///
	 scheme(s1mono) rows(2) saving(m4, replace) 
	
	
	
*Model 5
coefplot ///
	(matrix(m5_bl_sel[.,1]), se(m5_mid4_sel[.,2]) m(smcircle) mlw(medium) msize(medlarge) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m5_gml_sel[.,1]), se(m5_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m5_mid4_sel[.,1]), se(m5_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.2)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(col(3) ring()  region(lcolor(none)) size(2)) ///
	xtitle("Coefficient Value" , size(small)) ///
	ytitle(" ", size(small) margin(0 15 0 0) ) ///
	subtitle("Onset" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(-1 -.5 0 .5 1 , labsize(small) angle(0)  ) ///
	ylabel(none ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(20) xsize() fysize() ///
	graphregion(margin(l=0 r=18)) ///
	levels(95 90) saving(m5_sel, replace) 
coefplot ///
	(matrix(m5_bl_out[.,1]), se(m5_bl_out[.,2]) m(smcircle) mlw(medium) msize(medium) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m5_gml_out[.,1]), se(m5_gml_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m5_mid4_out[.,1]), se(m5_mid4_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.2)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(off) ///
	xtitle("") ///
	ytitle(" ", size(small) margin(0 15 0 0)) ///
	title("Model 5"  , size(medsmall)) ///
	subtitle("MID with over 250 Fatalities" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noline noextend titlegap(3) range(-1.25,1.25)) ///
	xlabel(none  ) ///
	ylabel(none) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(5) xsize() fysize(80) ///
	graphregion(margin(l=0 r=18 b=-8)) ///
	levels(95 90) saving(m5_out, replace)
graph combine m5_out.gph m5_sel.gph,  xsize(3) ysize(6) iscale(1.2) ///
	 scheme(s1mono) rows(2) saving(m5, replace) 

*Model 6
coefplot ///
	(matrix(m6_bl_sel[.,1]), se(m6_mid4_sel[.,2]) m(smcircle) mlw(medium) msize(medlarge) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m6_gml_sel[.,1]), se(m6_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m6_mid4_sel[.,1]), se(m6_mid4_sel[.,2])  m(smcircle) mlw(medium) msize(medlarge) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.2)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(col(3) ring()  region(lcolor(none)) size(2)) ///
	xtitle("Coefficient Value" , size(small)) ///
	ytitle(" ", size(small) ) ///
	subtitle("Onset" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noextend titlegap(3) range(-1,2)) ///
	xlabel(-1 -.5 0 .5 1 1.5 2 , labsize(small) angle(0)  ) ///
	ylabel(none ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(20) xsize() fysize() ///
	graphregion(margin(l=9 r=23.5)) ///
	levels(95 90) saving(m6_sel, replace) 
coefplot ///
	(matrix(m6_bl_out[.,1]), se(m6_bl_out[.,2]) m(smcircle) mlw(medium) msize(medium) mfcolor("253 174 97*.4") mlc("253 174 97")  ciop(lcolor("253 174 97" "253 174 97") lwidth(medthick thick) ) label(MID 3.0)) ///
	(matrix(m6_gml_out[.,1]), se(m6_gml_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("215 25 28*.4")  mlc("215 25 28") ciop(lcolor("215 25 28" "215 25 28") lwidth(medthick thick) ) label(GML Data)) ///
	(matrix(m6_mid4_out[.,1]), se(m6_mid4_out[.,2])  m(smcircle) mlw(medium) msize(medium) mfcolor("31 120 180*.4")  mlc("31 120 180") ciop(lcolor("31 120 180" "31 120 180") lwidth(medthick thick) ) label(MID 4.2)) ///
	,  scheme(s1mono)  ///
	coeflabels()   ///
	legend(off) ///
	xtitle("") ///
	ytitle(" ", size(small) margin(0 0 0 0)) ///
	title("Model 6"  , size(medsmall)) ///
	subtitle("War" , size(small)) ///
	ysca(noline titlegap(0)) ///
	xsca( noline noextend titlegap(3) range(-1,2)) ///
	xlabel(none  ) ///
	ylabel(none ) ///
	plotregion(color(gs15) style(none)) ///
	xline(0, lcolor(gray) lpattern(dash)) ///
	ysize(5) xsize() fysize(80) ///
	graphregion(margin(l=9 r=23.5 b=-8)) ///
	levels(95 90) saving(m6_out, replace)
graph combine m6_out.gph m6_sel.gph,  xsize(3) ysize(6) iscale(1.2) ///
	 scheme(s1mono) rows(2) saving(m6, replace) 
	
